function eigen_decomposition()
% EIGEN_DECOMPOSITION 特征分解详细演示
% 
% 展示特征分解的理论、计算和应用

fprintf('=== 特征分解详细演示 ===\n\n');

%% 1. 特征分解基础
fprintf('1. 特征分解基础\n');
fprintf('目标：理解特征分解的基本概念\n\n');

% 简单的2x2矩阵
A1 = [3, 1; 1, 3];
fprintf('矩阵 A1 =\n');
disp(A1);

% 计算特征值和特征向量
[V1, D1] = eig(A1);
eigenvals1 = diag(D1);
fprintf('特征值: [%.4f, %.4f]\n', eigenvals1);
fprintf('特征向量矩阵 V =\n');
disp(V1);

% 验证特征分解
A1_reconstructed = V1 * D1 * V1^(-1);
fprintf('重构误差: %.2e\n', norm(A1 - A1_reconstructed, 'fro'));

% 验证特征值方程 Av = λv
for i = 1:size(A1, 1)
    v = V1(:, i);
    lambda = eigenvals1(i);
    residual = norm(A1 * v - lambda * v);
    fprintf('特征值 %d 的残差: %.2e\n', i, residual);
end

%% 2. 对称矩阵的特征分解
fprintf('\n2. 对称矩阵的特征分解\n');
fprintf('目标：展示对称矩阵特征分解的特殊性质\n\n');

% 构造对称矩阵
A2 = [4, 2, 1; 2, 5, 3; 1, 3, 6];
fprintf('对称矩阵 A2 =\n');
disp(A2);

[V2, D2] = eig(A2);
eigenvals2 = diag(D2);

% 验证实特征值
fprintf('特征值（应该都是实数）:\n');
for i = 1:length(eigenvals2)
    fprintf('λ_%d = %.6f + %.6fi\n', i, real(eigenvals2(i)), imag(eigenvals2(i)));
end

% 验证正交性
orthogonality_error = norm(V2' * V2 - eye(size(V2)), 'fro');
fprintf('特征向量正交性误差: %.2e\n', orthogonality_error);

% 谱分解形式：A = λ₁v₁v₁ᵀ + λ₂v₂v₂ᵀ + λ₃v₃v₃ᵀ
A2_spectral = zeros(size(A2));
for i = 1:length(eigenvals2)
    v = V2(:, i);
    lambda = eigenvals2(i);
    A2_spectral = A2_spectral + lambda * (v * v');
end

fprintf('谱分解重构误差: %.2e\n', norm(A2 - A2_spectral, 'fro'));

%% 3. 非对称矩阵的特征分解
fprintf('\n3. 非对称矩阵的特征分解\n');
fprintf('目标：处理可能有复特征值的情况\n\n');

% 构造非对称矩阵
A3 = [1, 2, 0; 0, 1, 2; 2, 0, 1];
fprintf('非对称矩阵 A3 =\n');
disp(A3);

[V3, D3] = eig(A3);
eigenvals3 = diag(D3);

fprintf('特征值:\n');
for i = 1:length(eigenvals3)
    if abs(imag(eigenvals3(i))) < 1e-10
        fprintf('λ_%d = %.6f (实数)\n', i, real(eigenvals3(i)));
    else
        fprintf('λ_%d = %.6f + %.6fi (复数)\n', i, real(eigenvals3(i)), imag(eigenvals3(i)));
    end
end

% 验证特征分解
if all(abs(imag(eigenvals3)) < 1e-10)
    A3_reconstructed = real(V3 * D3 * V3^(-1));
    fprintf('重构误差: %.2e\n', norm(A3 - A3_reconstructed, 'fro'));
else
    fprintf('注意：存在复特征值，特征向量也可能是复数\n');
end

%% 4. 缺陷矩阵（不可对角化矩阵）
fprintf('\n4. 缺陷矩阵示例\n');
fprintf('目标：展示不可对角化矩阵的特征\n\n');

% Jordan块矩阵
A4 = [2, 1, 0; 0, 2, 1; 0, 0, 2];
fprintf('Jordan块矩阵 A4 =\n');
disp(A4);

[V4, D4] = eig(A4);
eigenvals4 = diag(D4);

fprintf('特征值: [%.6f, %.6f, %.6f]\n', eigenvals4);
fprintf('注意：特征值重数为3，但几何重数小于代数重数\n');

% 检查特征向量的线性无关性
rank_V4 = rank(V4);
fprintf('特征向量矩阵的秩: %d (应该小于 %d)\n', rank_V4, size(A4, 1));

% 计算条件数
cond_V4 = cond(V4);
fprintf('特征向量矩阵条件数: %.2e (很大表示接近奇异)\n', cond_V4);

%% 5. 特征值的扰动分析
fprintf('\n5. 特征值的扰动分析\n');
fprintf('目标：分析矩阵扰动对特征值的影响\n\n');

% 良态矩阵
A5_good = diag([1, 2, 3, 4]);
fprintf('良态对角矩阵:\n');
disp(A5_good);

% 病态矩阵
epsilon = 1e-6;
A5_ill = [1, 1, 1, 1; 
          0, epsilon, 1, 1; 
          0, 0, epsilon^2, 1; 
          0, 0, 0, epsilon^3];
fprintf('病态上三角矩阵 (ε = %.0e):\n', epsilon);
disp(A5_ill);

% 计算条件数
cond_good = cond(A5_good);
cond_ill = cond(A5_ill);
fprintf('良态矩阵条件数: %.2f\n', cond_good);
fprintf('病态矩阵条件数: %.2e\n', cond_ill);

% 添加小扰动
perturbation_level = 1e-10;
E = perturbation_level * randn(size(A5_good));

eigenvals_good_orig = eig(A5_good);
eigenvals_good_pert = eig(A5_good + E);
eigenvals_ill_orig = eig(A5_ill);
eigenvals_ill_pert = eig(A5_ill + E);

change_good = norm(sort(eigenvals_good_pert) - sort(eigenvals_good_orig));
change_ill = norm(sort(eigenvals_ill_pert) - sort(eigenvals_ill_orig));

fprintf('扰动水平: %.0e\n', perturbation_level);
fprintf('良态矩阵特征值变化: %.2e\n', change_good);
fprintf('病态矩阵特征值变化: %.2e\n', change_ill);
fprintf('敏感性比: %.1f\n', change_ill / change_good);

%% 6. 广义特征值问题
fprintf('\n6. 广义特征值问题\n');
fprintf('目标：求解 Ax = λBx 形式的广义特征值问题\n\n');

% 定义矩阵A和B
A6 = [2, 1; 1, 3];
B6 = [1, 0.5; 0.5, 2];

fprintf('矩阵 A =\n'); disp(A6);
fprintf('矩阵 B =\n'); disp(B6);

% 广义特征值分解
[V6, D6] = eig(A6, B6);
eigenvals6 = diag(D6);

fprintf('广义特征值: [%.6f, %.6f]\n', eigenvals6);

% 验证广义特征值方程 Av = λBv
for i = 1:length(eigenvals6)
    v = V6(:, i);
    lambda = eigenvals6(i);
    residual = norm(A6 * v - lambda * B6 * v);
    fprintf('广义特征值 %d 的残差: %.2e\n', i, residual);
end

% 与标准特征值问题的关系
if cond(B6) < 1e12  % B可逆
    A6_standard = B6 \ A6;
    eigenvals6_standard = eig(A6_standard);
    fprintf('通过B⁻¹A得到的特征值: [%.6f, %.6f]\n', eigenvals6_standard);
    fprintf('与广义特征值的差异: %.2e\n', ...
            norm(sort(eigenvals6) - sort(eigenvals6_standard)));
end

%% 7. 特征分解的应用：矩阵函数
fprintf('\n7. 特征分解的应用：矩阵函数\n');
fprintf('目标：使用特征分解计算矩阵函数\n\n');

A7 = [1, 0.5; 0.5, 2];
fprintf('矩阵 A7 =\n'); disp(A7);

[V7, D7] = eig(A7);
eigenvals7 = diag(D7);

% 计算矩阵指数 exp(A)
exp_eigenvals = exp(eigenvals7);
A7_exp = V7 * diag(exp_eigenvals) * V7^(-1);

fprintf('矩阵指数 exp(A7) =\n');
disp(A7_exp);

% 验证：与MATLAB内置函数比较
A7_exp_matlab = expm(A7);
fprintf('与MATLAB expm的差异: %.2e\n', norm(A7_exp - A7_exp_matlab, 'fro'));

% 计算矩阵平方根 sqrt(A)
if all(eigenvals7 > 0)  % 正定矩阵
    sqrt_eigenvals = sqrt(eigenvals7);
    A7_sqrt = V7 * diag(sqrt_eigenvals) * V7^(-1);
    
    fprintf('矩阵平方根 sqrt(A7) =\n');
    disp(A7_sqrt);
    
    % 验证：sqrt(A) * sqrt(A) = A
    verification = A7_sqrt * A7_sqrt;
    fprintf('验证 sqrt(A)² = A，误差: %.2e\n', norm(verification - A7, 'fro'));
end

% 计算矩阵对数 log(A)
if all(eigenvals7 > 0)
    log_eigenvals = log(eigenvals7);
    A7_log = V7 * diag(log_eigenvals) * V7^(-1);
    
    fprintf('矩阵对数 log(A7) =\n');
    disp(A7_log);
    
    % 验证：exp(log(A)) = A
    verification_log = V7 * diag(exp(log_eigenvals)) * V7^(-1);
    fprintf('验证 exp(log(A)) = A，误差: %.2e\n', norm(verification_log - A7, 'fro'));
end

%% 8. 特征分解的数值稳定性
fprintf('\n8. 特征分解的数值稳定性\n');
fprintf('目标：分析特征分解算法的数值特性\n\n');

% 测试不同条件数的矩阵
condition_numbers = [1e2, 1e6, 1e10, 1e14];
fprintf('条件数        特征值精度    特征向量精度\n');
fprintf('--------      ----------    ------------\n');

for cond_target = condition_numbers
    % 构造指定条件数的矩阵
    n = 4;
    eigenvals_target = [1, 1/sqrt(cond_target), 1/cond_target^(1/3), 1/cond_target];
    V_random = orth(randn(n, n));
    A_test = V_random * diag(eigenvals_target) * V_random';
    
    % 计算特征分解
    [V_computed, D_computed] = eig(A_test);
    eigenvals_computed = diag(D_computed);
    
    % 评估精度
    eigenval_error = norm(sort(eigenvals_computed) - sort(eigenvals_target)) / ...
                     norm(eigenvals_target);
    
    % 特征向量精度（通过残差评估）
    eigenvec_residuals = zeros(n, 1);
    for i = 1:n
        v = V_computed(:, i);
        lambda = eigenvals_computed(i);
        eigenvec_residuals(i) = norm(A_test * v - lambda * v) / (norm(A_test) * norm(v));
    end
    eigenvec_error = max(eigenvec_residuals);
    
    fprintf('%.0e      %.2e        %.2e\n', cond_target, eigenval_error, eigenvec_error);
end

%% 总结
fprintf('\n=== 特征分解总结 ===\n');
fprintf('1. 对称矩阵有实特征值和正交特征向量\n');
fprintf('2. 非对称矩阵可能有复特征值\n');
fprintf('3. 缺陷矩阵不能对角化\n');
fprintf('4. 特征值对矩阵扰动的敏感性与条件数相关\n');
fprintf('5. 广义特征值问题扩展了标准特征值问题\n');
fprintf('6. 特征分解可用于计算矩阵函数\n');
fprintf('7. 数值稳定性取决于矩阵的条件数\n\n');

end